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ABSTRACT 

We use SPH simulations to investigate the gravitational fragmentation of expanding 
shells through the linear and non-linear regimes. The results are analysed using spher- 
ical harmonic decomposition to capture the initiation of structure during the linear 
regime; the potential-based method of Smith et al. (2009) to follow the development 
of clumps in the mildly non-linear regime; and sink particles to capture the properties 
of the final bound objects during the highly non-linear regime. In the early, mildly 
non-linear phase of fragmentation, we find that the clump mass function still agrees 
quite well with the mass function predicted by the analytic model. However, the sink 
mass function is quite different, in the sense of being skewed towards high-mass ob- 
jects. This is because, once the growth of a condensation becomes non-linear, it tends 
to be growing non-competitively from its own essentially separate reservoir; we call 
this Oligarchic Accretion. 
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1 INTRODUCTION 



are ubiquitous 
on a wide 



Bubbles and the shells they sweep up 
features of the interstellar medium (ISM), 
range of mass an d size scales. IChurchwell et al.l (|2006l ) and 

IChurchwell et al.l (|2007l ) observed over 600 infrared shells 
with Spitzer on size scales of 0.1 to l.Opc. iKonvves et al.l 
■ (|2007T ) catalogued 462 lar ger-scale loops and arcs o bserved 
in the far infrared and lEhlerova fc PalousI (|2005f l found 
~ 1000 shell-like structures in the Leiden-Dwingcloo HI sur- 
vey on scales of tens to hundreds of parsecs. These shells are 
of intrinsic interest, since they form an important part of the 
structure of the ISM, but they are also interesting from the 
point of view of star formation, since it has long been sus- 
pected that dense shells are sites of triggered star formation; 
hence they form a cruc ial component in the self- propagating 
star formation model (|Elmegreen fe Ladalll977t ). 

There exists a considerable body of theoretical work 
on the fragmentation of expanding shel ls, usually i n the 
context of the thin-sh ell model fe.g lElmegreenl Il994l : 

IWhitworth et all Il994al lbl; IWtinsch fc PalousI 120011 ). The 
thin-shell model treats the expanding shell as infinitesimally 
thin and considers the linear behaviour of perturbations in 
the shell surface-density, treating the fragmentation as a lo- 
cally two-dimensional process. However, there is no clear 
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consens us on the mass function of stars produced by this 
process. IWhitworth et all <|l994bf ) assume that the most un- 
stable mode in the shell can be used to characterize the mass 
function, concluding that shell fragmentation should prefer- 
entially produce high-mass stars. If this is correct, it strongly 
supports the self-propagating star formation model, since 
each shell produces more than enough O -type stars to trig- 
ger th e next stellar generation. However, IWiinsch fc PalousI 
(2001) perform a more detailed analysis under the assump- 
tion that all unstable wavelengths are represented in the 
mass function. They derive a fragment mass function from 
the dispersion relation for the growth rate of perturbations 
of different angular wavenumbers, and conclude that shell 
fragmentation should produce a power-law mass function 
not very different from the canonical Salpeter function. 

Unambiguous observatio nal evidence for sh e ll frag men- 
tation is difficult to obtain. IChurchwell et alj (2007) find 
that ~ 18% of their sample of 269 bubbles in the inner 
Milky Way show evidence of triggered star formation, but 
attribute this to the shells over-running pre-existing density 
anomalies, rather than to fragmentation of the shells them- 
sleves. O bservations of s hells d riven by HII regions by , for 
example, Zavagno et alj (|2006l ) , IDeharveng et all |2006n and 
IDeharveng et al.l ~ 20081 ) often detect populations of young 
stars in the peripheries of bubbles and the stars are of- 
ten m assive. Herschel observations of Sh2-1 04 (|Rod6n et al.l 
|2010| ) and RCW120 (|Zavagno et all |2010| ) seem to imply 
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that shell fragmentation produces s mall numbers o f large 
fragments. However, observations by iDawson et al l (|2008l ) 
of the mass function of CO clumps in the Carina Flare su- 
pershell reveal a power-law mass spectrum. It is thus not 
clear from either an observational or a theoretical point of 
view how self-gravitating shells should fragment. 

The present paper is one of a series aimed at resolving 
these issues by performing detailed nu merical simulation s 
of fragmenting shells. In the first paper ijDale et al l l|2009h : 
hereafter Paper I) we studied the fragmentation of expand- 
ing shells in the linear regime using a Smoothed Particle 
Hydrodynamics (SPH) code and an Adaptive Mesh Refine- 
ment (AMR) code. We studied the growth of perturbations 
in the shell surface density and compared the results of the 
two codes with each other and with the thin-shell model. 
We found that the agreement between the two radically dif- 
ferent numerical schemes was excellent. However, we also 
found that the boundary conditions applied to the shell were 
of crucial importance in determining the mass-spectrum of 
perturbations. A shell expanding in a vacuum gets thicker 
as it expands, and this thickening suppresses fragmentation 
at short wavelengths, relative to the predictions of the thin 
shell model. Conversely, if the thickening of the shell is pre- 
vented by a confining pressure, the shorter wavelengths be- 
come unstable too and, for a certain finite value of the pres- 
sure, fragmentation proceeds in good agreement with the 
thin shell model. 

The second paper in this series, IWunsch et all (|20ld ) 
and hereafter Paper II, presents an alternative to the thin- 
shell model which allows the boundary conditions on the 
shell inner and outer surfaces to be included and thus the 
effects of external pressure to be treated analytically. We 
find that pressures higher than that required to keep the 
shell thickness constant extend the range of unstable wave- 
lengths towards lower values. This should in principle result 
in a fragment mass function with a steeper slope and we 
investigate this possibility below. 

In this paper, we extend some of the simulations of Pa- 
per I into the non-linear regime and allow the fragmentation 
of the shells to proceed to the formation of bound com- 
pact objects. These can be thought of either as protostars 
or protoclusters, depending on their masses, allowing us (i) 
to derive a mass function for the stellar/cluster populations 
formed by fragmenting shells; and (ii) to evaluate the ability 
of the thin-shell model to predict these mass functions. In 
Section 2, we discuss the thin-shell model and the means 
by which mass functions may be inferred from it. In Section 
3, we explain our numerical methods and in Section 4, we 
briefly outline the simulations conducted. Our results are 
presented in Section 5 and discussed in Section 6. Finally, 
we summarize our conclusions in Section 7. 



of spherical harmonic wavenumber I, 



2 MASS FUNCTIONS FROM THE THIN 
SHELL MODEL 

For a self-gravitating shell with instantaneous radius, ex- 
pansion velocity and mean surface density given by R, V and 
So, the infinitesimally thin shell model yields an expression 
for the growth rate u) of perturbations in the surface density 
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where c s is the sound speed inside the shell, and A and B are 
constants whose values depend on whether the shell sweeps 
up mass as it expands. In the case of a shell expanding into a 



vacuum, A 



B 



In the shells studied in Paper I and 



in this paper, the expansion velocity V of the shell is small 
compared with c s and so the first two terms in Equation [T] 

are small; 

In IWunsch et~ai1 (|2010l ). we derived an alternative ex- 
pression for the dispersion relation which drops the assump- 
tion that the shell is thin and takes into account the effect 
of the external pressure, Pext 
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where zq and ro are respectively the vertical half-thickness 
and azimuthal radius of a pe rturbation and eisa normalisa- 
tion parameter, described in IWunsch et all (|2010h . In order 
to be useful in the study of star formation, these dispersion 
relations must be converted into fragment mass functions. 

Be ginning from the thin-shell model, IWhitworth et all 
(|l994bh neglect the shell expansion terms and consider only 
the interplay between gravity and gas pressure within the 
shell and the masses of fragments formed. They consider 
perturbations of radius r, compute the acceleration g of ma- 
terial inside them, and write the timescale on which such a 
fragment condenses out of the shell as 
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IWhitworth et all l|l994bh then find the value of r that min- 
imises Equation [3] and take the mass corresponding to 
this fragm e nt siz e as representative of the mass function. 
iDale et all (|2007l ) simulate the early stages of the fragmen- 
tation of a shell driven by an expanding HII region a nd 
find reasonable agreement with I Whitworth et al.l (|l994bl ) in 
terms of the time and shell radius at which fragmentation 
begi ns and the mean masses o f fragments formed. 

IWunsch fc Palous! |200ll ) introduce a more sophisti- 
cated analysis in an attempt to derive the mass spectrum 
of fragments, but their derivation is not correct, since there 
is an error in their Equation 52. Although we import some 
of their results, we provide here a different semi-analytic 
derivation. 

We seek an expression for the number of fragments dN 
existing at a time t in the range of masses [m{,mi + Ami] 
and we make the assumption that the right-hand side of this 
expression can be separated into two parts; the fragmenta- 
tion integral (now written as a function of mass) 7f(mf,i)), 
which expresses how much perturbations of a given mass 
have grown at time t, and a geometrical function (J(mf)dmf 
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Figure 1. Result of an experiment of covering a square with 
circles uniformly drawn in mass from [1.26 x 1(T 5 ,7.9 x 1CT 1 ] 
up to a covering factor of 0.8. 




Figure 2. Mass functions of circles, with covering factors / of 
0.70 (plus signs), 0.75 (crosses) and 0.80 (triangles). The straight 
line has a logarithmic slope of -2.0. 



which encodes how many fragments of a given mass will fit 
on the shell. We may then write 



dN(m{,t) = Ii(mt,t)m { 20 Arm. 



(6) 



dN(nif,t) = It(mt, t)Q(rrif) Amj. 

The fragmentation integral is simply defined as 

If(m,{,t) = / u(mf , t')dt', 



(4) 



(5) 



where w(mf, t') is the instantaneous growth rate of the mode 
corresponding to mt and to is the time at which fragmenta- 
tion begins. 

The derivation of the function Q(m,j ) dmf is the point at 
which we depart from the methodology of lWiinsch fc Palousl 
l|200ll ). It is not clear to us that an analytic expression for 
G(m{)&mi may be found since, on the reasonable assump- 
tion that fragments may not overlap each other, the number 
of objects of any given wavenumber that may be accommo- 
dated by the shell depends on how much space has already 
been consumed by other fragments. We therefore adopted 
a Monte Carlo approach in an attempt to derive a semi- 
analytic approximation to 5(mf)dmf. 

We assume that fragments are circles. We took a unit 
square (the shape of the area to be filled is largely irrelevant, 
since circles do not tesselate) and populated it with circles 



chosen uniformly in [m m i n 



mf being the mass of a 



circle and proportional to its area. We forbade overlapping 
and continued until the square was covered up to a factor 
/ by circles. We then constructed histograms of the popu- 
lation of circles. In Figure [T] we show the result of one such 
experiment in which / = 0.8 and m m i n , m max are 1.26x 10~ 5 
and 7.9xl0 _1 respectively. In Figure [21 we show the power 
laws resulting from hundreds of such experiments in which 
we varied /. Varying the covering factor affects the low-mass 
end of the mass spectrum, which turns over at small masses 
due to underrepresentation of small circles. Achieving higher 
covering factors requires more small circles and pushes the 
turnover towards smaller masses. At the high-mass end of 
the mass function, the mass function becomes noisy. How- 
ever, in between these limits, the slope of the function is 
very robust. We find that the population of circles can be 
well approximated by diV oc mf 2 dm{. 

Armed with this result we may write 



An analytic expression for the fragmentation integral 
itself can be derived for the infinitesimally thin shell and is 
given in Appendix A in a dimensionless form that can be ap- 
plied to any thin shell. However, it is somewhat cumbersome 
and we compute fragmentation integrals from numerically 
integrating Equation [S] using analytic wavenumber growth 
rates for the thin-shell and PAGI models and compare the 
result to that computed from our SPH simulations. 



3 NUMERICAL METHODS 

We make use of the same SPH co de used in Paper I, a variant 
of that described i nl.Be.nq (119901) but more recently updated 
and described in iBate et ah \ 19951 ). The fluid equations 



are solved using the SPH technique implementin g the stan- 
dard a rtificial viscosity prescription described in MonagharJ 
(1992), with a = 1, = 2. The self gravity of the gas is 
included using a binary tree. Crucially for these calculation, 
the code allows very high density region s to be replaced b y 
point-mass sink particles, as described in lBate et al.l l|l995h . 

Once the density of a particle exceeds a threshold, set 
to 10~ 19 g cm -3 in these simulations, it and its ~ 50 neigh- 
bours are considered candidates for sink creation. In order 
for a dense particle and its neighbours to be replaced by a 
sink, the group of particles must (a) constitute more than a 
thermal Jeans mass; (b) be contracting; (c) be bound. Gas 
particles may subsequently be accreted by sink particles if 
they pass within the sink particle's accretion radius and (a) 
are bound to the sink; (b) are more bound to that sink than 
to any other sink; (c) do not have enough angular momen- 
tum to achieve a circular orbit around the sink. The mass 
resolutions in the medium- and high-pressure simulations 
presented here are 1.67M© and O.4M0 respectively. We use 
a sink accretion radius of 0.05pc which ensures that the ac- 
cretion flows at the accretion radius are always supersonic, 
so that accretion of gas particles does not produce spurious 
drops in pressure. We repeated part of the medium-pressure 
calculation using an accretion radius of 0.025pc, but found 
that this made no difference to the results. The code also 
allows for corrections to be made to the pressure and den- 
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sity gradients in the neighbourhood of sink particles to com- 
pensate for the fact that there are few or no gas particles 
inside the accretion radius, but we found that they had no 
effect on our results. In these calculations, the gas flows are 
non-turbulent and thus rather quiet and smooth, and hence 
comparatively easy to model. 

To achieve greater accuracy, the sink particles are not 
included in the binary tree and all gravitational forces in- 
volving them are computed by direct summation. To prevent 
sink particles in very close proximity to each other acquiring 
very short timesteps and stalling the simulation, we allowed 
point masses to merge if they passed within 0.05pc of each 
other and were mutually bound. During the course of the 
simulations, there are a few tens of such mergers, too few to 
strongly influence the mass functions produced. 

We use the same spherical harmonic decomposition 
technique discussed in Paper I to analyse the growth of per- 
turbations in the shell surface density, which we compare to 
the thin shell model. The thin-shell model uses linearized 
equations and is therefore only guaranteed to be valid in 
the linear phase of fragmentation, where surface density per- 
turbations are less than or comparable to the mean surface 
density. In this work we are interested in the transition of the 
shell into the non-linear regime and the contraction of such 
perturbations into bound objects. Apart from this phase of 
fragmentation being non-linear, there also comes a point 
for each contracting fragment when its wavelength becomes 
comparable to the thickness of the shell and fragmentation 
transitions from a two-dimensional to a three-dimensional 
process. In addition, the wavelengths and wavenumbers of 
fragments changes as they collapse, and the surface density 
associated with the fragment grows, so that the association 
of a given wavenumber with a given mass may become in- 
valid. For these reasons, it is not clear that the thin-shell 
or PAGI analysis may be extended into this regime of the 
shell's evolution and we therefore turn to a numerical study. 

We use two different techniques to follow the fragmen- 
tation process beyond the linear phase. Once fragmentation 
has become strongly non-linear and bound objects are col- 
lapsing on th eir local fre efall times, we make use of SPH 
sink particles (jBate et al.lll995l ). However, in the regime be- 
tween the linear and strongly non-linear phases of the shell 
evolution, before any fragments have become gravitationally 
bound, we require a robust clump-finding method. 

Three dimensional clump-finding in complex density 
fields is a difficult task and several numerical methods 
have been developed to facilitate it, e.g. the CLUMPFIND 
l|Williams et al.lll994 ) algorithm. In a particle-based code 
such as SPH, progress can be made in identifying perturba- 
tions in the density field by identifying contiguous regions 
of fluid whose density exceeds some threshold. Constructing 
a mass function is problematic however, since the chosen 
density threshold effectively selects the masses of the frag- 
ments. This problem can be solved by insisting that frag- 
ments must be bound. By selecting the densest particles as 
seeds for clumps and continually adding adjacent particles 
until the total energy of the composite object becomes posi- 
tive, it is possible to identify a population o f bou nd objects. 

This method was used by iDale et all l|2007l) . but is of 
little use here, since it is of no help in the mildly non-linear 
regime. We find that a given object goes very rapidly from 
becoming bound to forming a sink particle, so that this 
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Figure 3. Comparison of the PAGI low-pressure (1 X 10 17 dyne 
cm , green), PAGI medium-pressure (1 X 10 -13 dyne cm , 
black), PAGI high-pressure (5 X 10 — 13 dyne cm -2 , magenta), 
thin-shell (blue) and numerical (red) medium— pressure fragmen- 
tation integrals computed up to a time of 8.80 Myr. 



method reveals little that cannot be learned from the sink 
particles themsel ves. We therefore m ake use of the clump- 
finding method o f lSmith et aD (2009) which identifies struc- 
ture using peaks in the gravitational potential field, instead 
of in density. There are two advantages to such a method. 
First, the gravitational potential distribution in the cloud is 
considerably smoother than the density distribution, since 
density fluctuations that do not carry sufficient mass can- 
not contribute significantly to the potential field. Second, the 
strength of the gravitational potential determines whether a 
clump will collapse and how mass will flow. The structures 
identified by this algorithm are called p-cores to distinguish 
them from conventional density cores. 

The potential-clumpfinding algorithm works in a simi- 
lar manner to CLUMPFIND and can be applied directly to 
the SPH data. The SPH particle with the deepest gravita- 
tional potential forms the head of a clump, then the particle 
with the next deepest potential is either assigned to the 
same clump if it is an SPH neighbour of the head particle, 
or forms a new clump if not, and so on. Clumps are defined 
down to either a minimum positive potential, or the low- 
est contour which they share with a neighbouring clump. 
Unlike the traditional CLUMPFIND algorithm, we use con- 
tour levels primarily to define the level at which potential 
clumps join, rather than to distinguish clumps from noise, 
and so our contours are numerous and finely spaced. This 
has the effect of subtracting the background potential. This 
is necessary as gravity is a long range force, and so is af- 
fected by both the mass inside the p-core and surrounding 
it; hence we must remove the background to obtain the net 
effect on the mass within. P-cores, therefore, represent the 
local maximum above the surrounding background. 
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Figure 4. Comparison of the PAGI low— pressure (1 X 10 — 17 
dyne cm -2 , green), PAGI medium-pressure (1 X 10~ 13 dyne 
cm -2 , black), PAGI high-pressure (5 X 10~ 13 dyne cm -2 , 
magenta), thin— shell medium-pressure (blue) and numerical 
medium-pressure (red) mass functions, derived from Equation 
[6] computed up to a time of 8.80 Myr. The dashed line has a 
logarithmic slope of -2.25. 
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Figure 5. Comparison of clump mass functions generated by 
imposing an artificial surface density profile computed from the 
fragmentation integral of a given SPH dump on a uniform shell 
(red) with that computed from analysing the same dump directly 
(blue). 

4 SIMULATIONS 

In Paper I, we presented the results of simulations of 
expanding momentum-driven shells with two different 
boundary conditions. The shells had a mass of 2 x 10 4 M Q , 
an initial radius of 10 pc and an initial velocity of 2.1 
kms - , such that they expand to a maximum radius of 



« 23pc before beginning to contract. In the simulation 
with free boundaries, we observed that the shell thickened 
considerably during its expansion and that this suppressed 
the growth of fragments at large wavenumbers. We repeated 
the simulation applying constant and equal pressures of 
1.0 x 10~ 13 dyne cm -2 to the inner and outer faces of 
the shell such that its thickness remained approximately 
constant and found that this resulted in much better 
agreement with the thin shell model. However, in Paper I, 
we followed the shell evolution only in the linear regime in 
which the thin-shell model is most likely to be applicable. 
In this paper, we allow the simulations of Paper I to run 
further, into the non-linear regime, using only the SPH 
simulations and not the AMR simulations, since both sink 
formation and clumpfinding are much more difficult in 
the latter. We find that the fragmentation process in the 
non-pressure confined simulation from Paper I proceeds 
too slowly to form a significant number of sink particles 
before the shell has contracted back to its original size and 
we do not examine it further. 

In Paper II, we considered the effects of higher pres- 
sures and found that they result in additional shorter 
wavelengths becoming unstable, which should produce 
more low-mass fragments. In this work, we concentrate 
on the pressure-confined model from Paper I (which we 
refer to as the medium-pressure model), and we also study 
the mass function produced by the same shell, but with a 
confining pressure of 5.0 x 10~ 13 dyne cm* 2 , which we also 
studied in Paper II and will refer to as the high-pressure 
model. 



5 RESULTS 

5.1 Spectral analysis vs. clumpfinding 

In Figure [3] we plot the fragmentation integrals computed 
for the medium-pressure run from the thin-shell model, 
for three different pressures using the PAGI model, and 
that derived from the spherical-harmonic analysis of the 
medium-pressure SPH simulation, all computed at a time of 
8.80Myr. The analytic fragmentation integrals are obtained 
by integrating Equation [S] using Equation Q] to compute 
uj in the thin-shell case, and Equation [2] to compute ui ( 
in the PAGI case. The numerical fragmentation integral 
clearly agrees much better with the relevant PAGI model 
than with the thin-shell model. Importing the results 
from Section 2, we then compute mass functions from the 
fragmentation integrals, shown in Figure [4] We again see 
that the mass function computed from the PAGI model is 
in much better agreement with that computed from the 
SPH medium-pressure calculation than is the thin-shell 
mass function. We find that, at high masses, the theoretical 
mass functions are power laws with slopes of -2.25, in good 
agreement with the canonical Salpeter slope of -2.35. 

Before proceeding with clumpfinding, we sought to es- 
tablish a connection between the two-dimensional spectral 
analysis technique and the intrinsically three-dimensional 
concept of a clump defined by its gravitational potential. 
We selected an output file from the early stages of the 
medium-pressure simulation and computed the correspond- 
ing fragmentation integral. We then used the integral to 
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Figure 6. Comparison of mean (solid line) and maximum per- 
turbed (dashed line) surface density in the medium— pressure cal- 
culation. 

construct a perturbed surface density profile which we 
imposed on a uniform shell by iteratively adjusting the 
SPH particle masses in a manner similar to that described 
in Paper I. Finally, we computed a clump mass function 
from the artificial shell and compared it to that computed 
from the original output file from the medium-pressure 
simulation. As shown in Figure [5l the mass functions 
generated in this way are in reasonable agreement, which 
implies that the spectral analysis technique is a legitimate 
way of describing shell fragmentation, at least in the early 
linear regime. 



5.2 Mass function of clumps 

In Paper I, we performed tests in which we inserted 
monochromatic density perturbations into the shell and 
compared their evolution with that predicted by the thin 
shell model. We found that the agreement was very good 
until the point when the perturbation in the surface density 
became equal to the mean surface density, after which the 
perturbation growth became non-linear, departing strongly 
from the model. In Figure [6l we plot the mean shell sur- 
face density (solid line) and the maximum perturbation in 
the surface density (defined as the maximum surface density 
minus the mean surface density, dashed line). We see that 
the time at which the surface density perturbation comes 
to exceed the mean surface density, and therefore the point 
at which perturbation growth is expected to become non- 
linear, occurs at ta 13Myr. 

In Figure [7] we plot the mass functions of the p-cores 
detected by the clumpfinding procedure at three epochs in 
the medium- and high-pressure calculations, and the slope 
derived from integration of the analytical fragmentation in- 
tegral. We find that, in the early stages of the shell's evo- 
lution (left panels) while the shell's evolution is still in the 
linear regime, the mass spectrum of fragments identified by 
the clump-finder follows has a form very similar to that pre- 



dicted by the PAGI model, although the p-core mass func- 
tion appears to be offset somewhat towards lower masses. 
There are two reasons for this. Firstly, material on the out- 
skirts of a given perturbation is rejected by the potential 
clump-finding algorithm as being part of the background, so 
that the algorithm underestimates the masses of fragments. 
More importantly, the association of perturbations of wave- 
length A with objects of mass m/ requires an assumption 
about the radius r of the object in terms of the wavelength. 
We adopt r — A/2, but this likely overestimates the masses 
of fragments as it includes material whose density is little 
different from that of the background shell. 

At later times (but before the formation of the first sink 
particles) the mass function of the p-cores becomes top- 
heavy (centre and right panels of Figure 0, with a deficit of 
low-mass objects, although the slope of the high-mass end 
of the mass function remains consistent with a power law. 

In Figure [5] we plot the numbers of p-cores found by 
the clumpfinder (before the formation of any sink particles) 
against the boundedness of the clumps _E ra t, defined as 



where is the kinetic energy calculated with respect to the 
centre of velocity of the p-core, -Etherm is its thermal energy 
and E p is the potential energy of the p-core calculated using 
the relative depth of the potential well once the background 
has been subtracted. Thus a p-core gravitational potential 
depth is expressed relative to the potential level in its vicin- 
ity. Values of J5 rat larger than unity indicate p-cores bound 
with respect to their environment, not merely when consid- 
ered in isolation. Figures [6] [7] and [8] show a clear sequence of 
events in which the evolution of the shell begins to depart 
from the predictions of the linear PAGI analysis around the 
time when the surface density perturbations come to exceed 
the mean shell surface density. The clump mass function 
then ceases to be a power law and the perturbations begin 
to contract and become gravitationally bound. 

A fundamental assumption of the thin-shell model is 
that all fragments evolve independently of one another in an 
otherwise unperturbed shell. We used the potential clump- 
finding code to verify this assumption. The clumpfinding 
code is able to follow individual p-cores to see how they gain 
or lose mass and whether they merge or exchange mass with 
each other as time progresses. To make use of this facility, 
we arbitrarily divide the first 15Myr of the medium-pressure 
shell's evolution into ten contiguous 1.5 Myr epochs. In Fig- 
ure [9] we plot the total number of p-cores extant in each 
epoch together with the numbers of clumps formed and de- 
stroyed in the previous epoch, against time. Initially, the 
number of p-cores destroyed in the previous epoch is al- 
most equal to the total number of p-cores, implying that 
the objects detected by the clumpfinder are transient enti- 
ties. However, as the simulation progresses, the total num- 
ber of p-cores rises and the number destroyed falls, so that 
the objects clearly become more longlived. Around 15Myr, 
the total number of p-cores levels off (and the numbers of 
new ones formed and of existing ones dissolving declines al- 
most to zero) implying that the shell evolves to contain a 
stable number of long-lived objects. In Figure [10] we plot 
the numbers of p-cores in each 1.5Myr epoch which have 
survived from the previous epoch, the numbers which have 
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Figure 7. Mass functions of clumps identified by the potential— based clump finder in the medium— pressure calculation (top row) and 
in the high-pressure calculation (bottom row) at three epochs (histograms). Left panels show the mass functions derived from the PAGI 
model (solid lines) and centre and right panels show power— laws with a logarithmic slope -2.25 (dashed line). 
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Figure 8. Boundedness (as defined by Equation of fragments detected by the potential clump-finding algorithm at three epochs in 
the medium-pressure calculation. 



retained more, or less, than half of the particles previously 
assigned to them and the numbers that have merged since 
the previous epoch. The number of p-cores surviving from 
the previous epoch rises rapidly as they become longlived, as 
opposed to ephemeral. We find that, in the vast majority of 
cases an SPH particle assigned to a given p-core remains as- 
signed to that p-core unless the p-core itself dissolves, and 
that mergers or fragmentation of cores are very rare. Our 
findings therefore support the assumption that surface den- 



sity perturbations do not interact with each other. We also 
note that the total number of p-cores levels off at a value of 
~ 200. Figure [3] shows that the most unstable wavenumber 
at the time when the first p-cores become bound is ~ 25. If 
the mass of a fragment is 7rA 2 M/16ttR 2 , the number of frag- 
ments with this wavenumber that the shell can support is 
4l 2 /ir 2 ~ 250, in reasonable agreement with the total num- 
ber of fragments detected by the clumpfinder. 
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Figure 11. Mass functions of sink particles at three epochs in the medium-pressure calculation (top row) and the high-pressure 
calculation (bottom row). 



5.3 Sink particles and oligarchic accretion 

Once the p-cores start to become bound, they collapse and 
form sink particles and the mass function is more easily fol- 
lowed by tracing the masses of the sinks themselves, shown 
in Figure 1111 We see that the sink-particle mass functions 
rapidly come to resemble the clump mass functions pre- 
sented in the previous section, since there are, by the epoch 
of sink-formation, very few clumps that are not bound and 
collapsing. The sink mass functions become ever more top- 
heavy, and even the high-mass end of the mass function loses 
its power-law slope. We note that the medium-pressure 
mass function peaks at a mass of ~ 8OM0 , corresponding to 
the mass associated with the most unstable wavenumber of 
I ~ 25 from Figure 

Implicit in the thin-shell model is the assumption that 
a single perturbation will produce a single collapsed object 
whose mass is proportional to the mass of the perturbation, 
which in turn is given by A 2 E/16, where A is the perturba- 
tion wavelength and E is the mean shell surface density at 
the time when the fragment collapses. In this picture, the 
accretion rate onto each collapsed object thus rises to a peak 
and falls back down to zero. On the contrary, as shown in 
Figure 1121 we find that, although the accretion rate onto 
a given sink is a maximum at the time of the sink parti- 
cle's formation and the accretion rate tails off, for most sink 
particles it never falls to zero and they continue accreting 
for the duration of the simulation so that, as shown in Fig- 
ure [131 sink masses are approximately proportional to their 



ages. We find that for most sink particles, the time-averaged 
accretion rate is a strong function of neither age nor mass. 
The uniformity in the accretion rates of the sink particles 
is a consequence of the fact that most of them form in the 
period of time when the shell is almost stationary at its max- 
imum radius, so that most sinks are born embedded in gas 
of the same density. Since the shell and the sinks are not 
in relative motion, the rate of accretion onto a given sink is 
determined only by the local freefall time, which is in turn 
the same throughout the shell at any given time. 

In Figure [TH we show the rate of sink particle formation 
and the cumulative number of sink particles as functions of 
time in the simulation. The number initially rises strongly 
before beginning to tail off around 23 Myr. We showed in 
Figure [9] that, as the first p-cores become bound, the to- 
tal number of p-cores stabilises. The final number of sink 
particles in Figure [TJ] is very similar to the final number of 
p-cores, and the tailoff in the rate of sink-formation appears 
to be due simply to the fact that there are only a finite num- 
ber of cores available from which to form sink particles. The 
top-heavy form of the mass function is then a consequence 
of the early pulse of sink-particle formation coupled with 
the fact that all sink particles share a common environment 
and thus accrete at comparable rates. We refer to this pro- 
cess as oligarchic accretion, since it is the sink particles that 
form first that accrete most of the shell's mass. 
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Figure 9. Total number of p— cores found by the clumpfinder 
in the medium-pressure run as a function of time (solid line), 
number of p— cores destroyed in the previous epoch (dashed line) 
and number of new p-cores formed in the previous epoch (dotted 
line). Epochs are 1.5Myr in duration.. 



6 DISCUSSION 

The purpose of this paper was to simulate the fragmentation 
of an expanding shell, following the evolution well into the 
non-linear regime to see if the predictions of the thin-shell 
and PAGI models (nominally only valid in the linear regime) 
can be extrapolated to accurately determine the mass func- 
tion of fragments formed. We find that the mass function 
of objects located by our clump-finding code agrees reason- 
ably well with the predictions of the thin-shell model in the 
early stages of fragmentation before significant numbers of 
fragments become bound. We also find that the assumption 
implicit in the model that fragments do not interact appears 
to be sound. However, once p-cores start to become bound, 
the numbers of p-cores detected levels off at a number con- 
sistent with the most unstable wavenumber at this epoch. 
Once the first cores become bound, they suppress the for- 
mation of new fragments. Once sink particles begin to from 
from the p-cores, the rate of sink formation is initially high, 
but tails off as the available cores are consumed. Sink par- 
ticles continue to accrete at roughly constant rates, since 
they all inhabit a similar environment, so that the initial 
pulse of sink-particle formation translates into a top-heavy 
mass function in which the mass corresponding to the most 
unstable wavelength at the time of the p-cores becoming 
bound is o ver-represented w i th resp ect to other modes. The 
analysis of lWhitworth et al.l (|l994bl ). in which the first un- 
stable wavenumber is taken as representative of the mass 
fu nction, describes fragmen tation better than the analysis 
of lWunsch fc Palousl l|200ll ) in which all modes contribute 
to the mass function. 

We refer to this process as oligarchic accretion, since 
the objects that form first win simply by virtue of being 
first. This mechanism i s very different from t he competitive 
accretion described by iBonnell et alj ( |2001 ;) , since in that 
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Figure 10. Number of p-cores surviving from previous epoch 
(solid line), number of p-cores retaining more than half of the 
particles assigned to them in the previous epoch (dashed line), 
number of p— cores retaining less than half of the particles assigned 
to them in the previous epoch (dotted line) and number of p-cores 
involved in mergers in the previous epoch (dash-dotted line), all 
plotted as functions of time, and all from the medium-pressure 
run. Epochs are 1.5Myr in duration. 
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Figure 12. Evolution of the mass of a randomly chosen 10% of 
the sinks formed in the simulation as a function of time in the 
medium-pressure calculation. 



model it is the different environments of the few stars in the 
dense gas at a cluster centre that enable them to accrete 
more than their siblings. In the simulations presented here, 
most of the sink particles form in gas of the same density, 
around the time when the shell expansion stalls. It is instead 
the time of formation and the non-uniform rate of formation 
of sink particles which is important in determining the final 
mass of a given object. Those objects forming first are able 
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Figure 13. Masses of sink particles as a function of their age as 
measured at a time of 31 Myr in the medium-pressure calculation. 



a 




22.5 25 
Time (Myr) 

Figure 14. Rate of formation (solid line) and total number 
(dashed line) of sink particles in the medium— pressure run plotted 
against time . 



to accrete more mass. The mass function becomes skewed 
because the rate of sink particle formation declines as the 
supply of p-cores forme d during the linear phase o f the shells 
evolution is consumed. iKlessen fc Burkertl (|200Gh observe a 
similar process in the evolution of a turbulent protocluster, 
where the mass function becomes skewed towards higher 
masses as the star formation efficiency becomes large and 
the gas reservoirs required to form new cores are depleted. 
The accretion process may also be aided by the fact that, at 
a time of «20 Myr, the shell begins to contract. However, 
the contraction is slow compared to the rate of sink particle 
formation and accretion, and the fragment mass function is 
already strongly top-heavy at this epoch, so contraction of 
the shell cannot be the primary driver of oligarchic accre- 
tion. 



It is possible that, in the case of a shell sweeping up mass 
as it expands, so that reserves of fresh gas in the shell are 
constantly replenished, p-cores would be able to form at all 
times, resulting in a mass function more closely resembling 
a power law. Howev er, such a shel l would experience the 
Vishniac instability (| Vishniad 1 19831 ) . probably altering the 
mass spectrum of fragments in the non-linear portion of the 
shell's evolution. This is important, since our results suggest 
that the most unstable mode at the time when bound ob- 
jects begin to form is over-represented in the mass function. 
To unequivocally demonstrate that accretion of fresh mate- 
rial would allow unabated p-core formation and produce a 
power-law mass function would require that the shell's gas 
reservoir be replenished in some way during sink formation. 
This is difficult in the case presented here, since most of the 
sinks form when the shell is almost stationary. The only way 
to replenish the gas would then be to actively feed matter 
into the shell, perhaps allowing it to infall from a reservoir 
just outside the shell's maximum radius, a highly artificial 
construction. We therefore defer answering this question to 
a later paper in which we analyse the effect of the Vishniac 
instability on a momentum-driven shell sweeping up an ex- 
ternal medium. 

The mass functions produced by our simulations bear 
little resemblance to any known stellar or cluster mass func- 
tion. As explained in Section 3, the mass resolution of our 
sink particles is good enough to resolve low-mass stars, but 
the linear resolution implied by the accretion radii (0.05pc) 
is rather large. Our sink particles should strictly be regarded 
as binaries or small multiple systems so that our mass func- 
tions, in the terminology of ?, are Multiple Star Mass Func- 
tions (MSMFs), not Single Star Mass Functions (SSMFs). 
This would be true to an even greater degree if our results 
were rescaled to higher masses. ? find, not surprisingly, that 
the SSMF contains fewer high-mass objects and more low- 
mass objects relative to the MSMF, and that the SSMF 
peaks at lower masses. The SSMF from our simulations 
would probably therefore be slightly less deficient in low- 
mass objects and contain fewer high-mass objects than the 
mass functions we present, but the effects of multiplicity will 
not quantitatively change our finding that the mass func- 
tion s become very t op-he avy due to oligarchic accretion. 

iDawson et al.l (|2008l ) observed the mass function of 
CO clumps in the Carina Flare supershell and observed 
that their mass spectrum follows a power-law. Although 
they caution against over-interpretation of their data, their 
clump mass measurements imply that the clumps are not 
gravitationally bound (in the sense that their inferred 
masses are below their virial masses). This result is con- 
sistent with our findings that the as-yet-unbound objects 
produced in the linear phase of shell fragmentation have 
a power law mass function (although note that the power 
law we derive i s con siderably steeper than that inferred by 
IDawson et al.l l|2008l )). Our findings suggest, however, that 
the mass function of these unbound fragments cannot nec- 
essarily be used to predict the mass function of stars or star 
clusters that will eventually form. There are as yet no obser- 
vations of the stellar or cluster mass functions generated by 
the fragmentation of expanding shells, owing to the extreme 
difficulty of making such a measurement. 

The initial conditions used in our simulations, that of 
a shell which is pressure-confined and yet does not sweep 
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up any mass are somewhat unusual and the peculiar mass 
function we obtain implies that most stars or clusters do 
not form in shells of this type. However, such shells must 
surely occur sometimes, since a shell may expand into the 
hot rarefied ISM and be pressure-confined while accreting 
negligible mass. Our results demonstrate that shells of this 
nature will produce anomalously high numbers of massive 
stars or clusters and may therefore lead to the propagation 
of star formation by triggering. 



7 CONCLUSIONS 

We find that, in the case of a momentum-driven shell 
pressure-confined internally and externally, the PAGI model 
faithfully reproduces the mass function of fragments pro- 
duced by the gravitational instability during the time for 
which perturbations in the shell surface density are rela- 
tively small and the behaviour is linear. 

In the case studied here of a shell which does not sweep 
up mass, we find that a form of accretion which we term 
oligarchic accretion operates within the shell whereby the 
first objects to form become the most massive because they 
accrete from the gas reservoir for longer. Coupled with the 
finite number of p-cores formed during the linear fragmenta- 
tion of the shell, this leads to a very top-heavy mass function 
in which the most unstable wavelength at the time when per- 
turbations first become bound is over-represented. A shell 
whose gas reservoir is replenished may form fragments con- 
tinuously, resulting in a mass function more closely resem- 
bling a power law. 

Since the mass function generated by our simulations 
is so unusual, we infer that shells of the type we have stud- 
ied cannot make a strong contribution to the global stellar 
population, although no direct measurements of the mass 
functions produced by expanding shells exist against which 
we could test our conclusions. 
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9 APPENDIX A 

The thin shell dispersion relation 
AV 
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may be recast in a dimensionless form with the help of an 
expansion law - a prescription for the variations in the shell 
radius R, expansion velocity V and surface density E with 
time. The shells described in this work are ballistic and of 
fixed mass M, expanding due to their inertia into a vacuum 
and decelerated by self-gravity alone. 

We define a dimensionless radius £ as 
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where K = MV 2 /2 and W = (GM 2 )/(2R) are the kinetic 
and potential energy respectively, and 
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is the maximum radius to which the shell expands (the initial 
total energy Ko + Wo at radius Ro is assumed to be negative 
so that the shell is bound). The free fall time obtained from 
Kepler's third law is 
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t S = JJ^ n (11) 
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and the most unstable dimensionless wavenumber at -R max 
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Inserting ((9} and (1121) into |[5]| and normalizing it by (|11[) we 
obtain 
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Using numerical factors A — 3/2 and B = 1/4 of a non- 
accreting shell we obtain the final dispersion relation 



w(Z,£, L 



.*E ( l_ Va r 3/ a . 



C (1 + 2/) -c 1 + 



In 



I 1/2 



(14) 



The degree of fragmentation at a given wavenumber I is 
described by the fragmentation integral 
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where equations © and JTT} have been used, we may write 
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Inserting (|18ft into (|15|l yields the fragmentation integral in 
a form 
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Using the Sage symbolic manipulation package (available at 
www.sagemath.org), the fragmentation integral becomes 
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